function [err, EPSCresult]=EPSC_Rect()
% Analyze the rectification of EPSPs
% Derived from EPSP_TC
% Expects data to be an EPSC, with alternating stim (stim on/stim off)
% over a range of different voltages. Subtracts the currents out and 
% plots voltage dependence of EPSC as well as subtracted traces.
%
% Paul B. Manis, Ph.D.
% pmanis@med.unc.edu
%
% Initial versions: 4/2002, 5/2002.
% 6/2/2002.
% 11/25/02
% Revised 8/18/04. P. Manis
% Interpolate IV's to common scale (that of the EPSC data) and 
% make calculations based on corrected traces.
% Added Woodhull fitting to match polyamine (or other ion) block.
%
persistent def_pA def_T


try
   sf = getmainselection;
   if(sf > 0) 
      pflag = getplotflag;
      QueMessage('EPSC_Rect analysis', 1); % clear the que
      % get some general parameters that we need to include for the
      % analysis
      if(isempty(def_pA))
          def_pA= 50;
          def_T = 22;
      end;
      def = {sprintf('%d', def_pA), sprintf('%d', def_T)};
      prompt = {'[polyamine] (uM): ', 'Temp (deg C): '};
     dlgTitle='Input rectification analysis';
     lineNo=1;
        answer=inputdlg(prompt,dlgTitle,lineNo,def);
    pA = number_arg(answer{1});
    T = number_arg(answer{2});
     def_pA = pA;
     def_T = T; % save in case we run again
     
      for i = 1:length(sf)
         EPSC_Rect2(sf(i), pflag, T, pA);
      end;
   end;
catch
   watchoff;
   QueMessage('Error in EPSC_Rect Analysis routine', 1);
end;

%-----------------------------------------------------------

function [err, EPSCresult]=EPSC_Rect2(sf, plot_flag,  T, pA)

global VOLTAGE CURRENT DFILE
global CONTROL
% first generate empty output arrays

err = 0;
EPSCresult=[];

[DFILE, err] = analysis_setup(DFILE, sf); % do default setup.
if(err ~= 0)
   return;
end;

%-----------------------------Begin analysis--------------------------------
QueMessage('EPSC_Rect - Starting');
dat = [];
time = [];

% abstract general information
protocol=deblank(lower(CONTROL(sf).protocol));
rate = DFILE.rate.*DFILE.nr_channel/1000;
[records,pts]=size(CURRENT);

[RL, err] = record_parse(CONTROL(sf).reclist);
frec=min(RL);
lrec=max(RL);

QueMessage('EPSC_Rect - time base');
% compute the time base for plotting (time is [rec, npts] array for EACH record)
time=make_time(DFILE);
tmax=max(max(time));
k = find(DFILE.rate < 1);
DFILE.rate(k) = 50;
RATES = (DFILE.rate .* DFILE.nr_channel)/ 1000; % array of sampling rates, convert to msec

% always compute the time of each stimulus, in seconds from the start of the data
QueMessage('EPSC_Rect - time base 2');
if(DFILE.mode >= 5)
   wz=DFILE.ztime;
   w=find(diff(wz) < -12*60*60); % correct for possibility that someone actually records the data across midnight rollover.... (yes, it happened. 5/16/01 with Huijie's data.)
   if(~isempty(w))
      wz(w+1:end)=wz(w+1:end)+24*60*60;
   end;
   zt = (wz-wz(1))/(60);
   cond_baseline = 5; % 5 min baseline
else
   zt = (DFILE.ztime-DFILE.ztime(1))/(60*1000);
   cond_baseline = 3;
end;
TM=zt;

QueMessage('EPSC_Rect - time base 3');
if(DFILE.mode >= 5)
   wz=DFILE.ztime;
   w=find(diff(wz) < -12*60*60); % correct for possibility that someone actually records the data across midnight rollover.... (yes, it happened. 5/16/01 with Huijie's data.)
   if(~isempty(w))
      wz(w+1:end)=wz(w+1:end)+24*60*60;
   end;
   ZT = (wz-wz(1))/(60);
   cond_baseline = 5; % 5 min baseline
elseif length(DFILE.ztime) == length(RL)
   ZT = (DFILE.ztime-DFILE.ztime(1))/(60*1000);
   cond_baseline = 3;
else
   ZT = ones(length(RL), 1)*DFILE.cycle;
   cond_baseline = 5;
end;

% access ZT with the record number (1..n) to get the corresponding ztime
% have to get the records first

% Now, get the times when the valves switched (if any...)
% and generate periods with valve 1, 2, 3 or 4.
QueMessage('EPSC_Rect - Valves...');
p=datac('getnote'); % read the current notefile information.
t_sw_valve=[];
n_valve=[];
TL=[];
VL=[];
if(~isempty(p) & length([p.proto]) > 0) % there should be some, but if not, don;'t do much
   % first set of arrays are immediate representations.
   sw_valve=[1 [p(find(diff([p.valve])~=0)+1).frec]];	% valve switch list (records)
   if(length(sw_valve) > 1)
      n_valve=[1 p(find(diff([p.valve])~=0)+1).valve]; % which valve...
      t_sw_valve=ZT(sw_valve); % don't forget offset from start of data..
      % now make long time arrays to match the other arrays.
      for i = frec:lrec
         TL(i)=ZT(i);
         for j=1:length(sw_valve)
            if(i >= sw_valve(j))
               VL(i)=n_valve(j);
            end;
         end;
      end;
   end;
end;

% Get analysis windows 
QueMessage('EPSC_Rect - analysis windows');
stim_list=CONTROL(sf).stim_time; % get the array
psp_time=CONTROL(sf).psp_time; % get psp definition array
if(ischar(psp_time))
	psp_time = number_arg(psp_time);
end;
if(ischar(stim_list))
	stim_list= number_arg(stim_list);
end;
%l = length(psp_time);
%if(l > 4)
%	psp_time = reshape(psp_time, l/2, 2)';
%end;
p2=size(psp_time);
npsc = length(stim_list); % get number of pscs to examine

trmp = floor(number_arg(stim_list(1))./RATES - 0.5); % for RMP/ihold determination
for i = 1:npsc % compute time windows for each psc in the train
   t0(i,:) = floor(number_arg(psp_time((i-1)*2+1))./RATES);
   t1(i,:) = floor(number_arg(psp_time((i-1)*2+2))./RATES);
end;
tep1 = floor(number_arg(stim_list(1)-1)./RATES);
tep2 = floor(number_arg(stim_list(1)+10)./RATES);

%-----------------------------Begin analysis in earnest --------------------------------
QueMessage('EPSC_Rect - Measuring Ih and rmp');

% measure holding current and "rmp"
% also measure Rin with s1 pulse after ho (ts2)
% and tau by fitting to s1 pulse.
for i = 1:records
	hold_cur(i) = mean(CURRENT(i,1:trmp(i))); % save in array for later usage
	RMP(i) = mean(VOLTAGE(i,1:trmp(i)));
end;
CONTROL(sf).iHold = mean(hold_cur);
CONTROL(sf).Rmp = mean(RMP);

% Now smooth the Current out a bit
QueMessage('EPSC_Rect - Smoothing');
for i = 1:records
	fsamp = 1000/RATES(i); % get sampling frequency
	fco = 10000;		% cutoff frequency in Hz
	wco = fco/(fsamp/2); % wco of 1 is for half of the sample rate, so set it like this...
	if(wco < 1) % if wco is > 1 then this is not a filter!
	   [b, a] = fir_win(8, wco); % fir type filter... seems to work best, with highest order min distortion of dv/dt...
	   ismo(i,:) = DigitalFilt(b, a, CURRENT(i,:)')'; % filter all the traces...
	else
	   ismo(i,:) = CURRENT(i,:);
	end
end

QueMessage('EPSC_Rect - Finding EPSC peak currents');
% first get baseline current
for r = 1:records
   ibase(r)=mean(CURRENT(r,1:trmp(r)));
end;

% Now find EPSC amplitudes for each entry
k = 1;
idiff=zeros(floor(records/2),size(ismo,2));
truers = 0.001*CONTROL(sf).access*CONTROL(sf).comp/100
for r = 2:2:records
   idiff(k,:) = ismo(r,:)-ismo(r-1,:); % take local difference
   ie(k,:) = ismo(r,:);
   ik(k,:) = ismo(r-1,:);
   vmrec(k) = mean(VOLTAGE(r,t0(1,r):t1(1,r))); % this is voltage, but not the true voltage
   % calculate the true voltages for both sets of traces
   ve(k,:) = VOLTAGE(r,:); %-ie(k,:)*truers;
   vk(k,:) = VOLTAGE(r-1,:); %-ik(k,:)*truers;
   k = k + 1;
end;
nstep = k-1;

[vm, ko] = sort(vmrec); % get the order of the voltages, since they may be randomized in the acquisition

vk=vk(ko,t0(1,r):t1(1,r));
ik=ik(ko,t0(1,r):t1(1,r));
ve=ve(ko,t0(1,r):t1(1,r));
ie=ie(ko,t0(1,r):t1(1,r)); % sorts all the arrays!
% now that we have corrected voltages, we can interpolate to get the
% current in the control traces at the same voltages as the stimulated
% traces
% 
QueMessage('EPSC_Rect - Interpolating Voltages');
    for i = 1:size(vk,2) % just do the windowed segment we need
        iki(:,i) = interp1(vk(:,i), ik(:,i), ve(:,i), 'spline'); % but do over all points in the IV at each time
    end;
% iki contains the interpolated currents for the control traces    
for k = 1:nstep
   i1 = max(ie(k,:)-iki(k,:)); % idiff(r,t0(1,r):t1(1,r)));
   i2 = min(ie(k,:)-iki(k,:)); %(idiff(r,t0(1,r):t1(1,r)));
   if(abs(i2) > i1) 
      i1 = i2;
   end;
   IEPSC(k) = mean(ie(k,:)-iki(k,:)); % i1;
   IPSCnc(k) = mean(ie(k,:) - ik(k,:)); % uncorrected IPSC
   vm(k) = mean(ve(k,:));
end;
QEPSC = 0; % not calulated yet.
%[vm,k]=sort(vm);
%IEPSC = IEPSC(k);
QueMessage('EPSC_Rect - Polyfit');
p = polyfit(vm, IEPSC, 4);
Vp = polyval(p, vm);


[vmin, kmin] = findnear(vm, -60);
imin = polyval([4*p(1) 3*p(2) 2*p(3) p(4)], vmin); % evaluate slope
[vmax, kmax] = findnear(vm, +40);
imax = polyval([4*p(1) 3*p(2) 2*p(3) p(4)], vmax); % evaluate slope

rratio = imax/imin; % this is the ratio of the SLOPES!, not the currents...
vrev = roots(p);
[w, krev] = min(abs(vrev)); % find the root closest to 0.
vrev = vrev(krev);

RTF=25.3*T/22; % scale by temperature that has been passed.
Vr = vrev; % reversal
% fit a woodhull equation to the data also
QueMessage('EPSC_Rect - Woodhull fit');
warning off MATLAB:divideByZero
warning off MATLAB:nearlySingularMatrix
ic{1} = [5 50 pA 1];
ic{2} = [1 100 200 1];
ic{3} = [10 25 100 1];

for i = 1:length(ic) % do different initial conditions
    [ex{i}, px{i}] = curve_fitting(vm, IEPSC, 'levenberg', 'cubic', 75, 4, [ic{i}], [1 1 1 1], [0 -500 -300 -50], [1000 500 300 50], RTF, Vr, 200);
end;
[minerr, ierr] = min([ex{:}]); % find the best one
e = ex{ierr};
p = px{ierr};

ywoodhull = -fit_func(p, vm, zeros(1, length(vm)), 75, RTF, Vr);
 % (vm-Vr).*g0./(1+p(1)./(p(2)*exp(-p(3)*vm/RTF)));
fprintf(1, 'Woodhull Parameters: [Blocker] = %8.3f  Kd(0) = %8.3f   z(1-delta) = %8.3f\n', p(2), p(3), p(4));
fprintf(1, 'Rectification Index (KD(0)/[pa]) = %8.3f  Vr = %8.2f\n', p(3)/pA, Vr);
fprintf(1, 'Fit error: %8.3f\n', e);

%fprintf(1, 'IEPSC(+40) = %6.2f  V(+40)=%6.2f %6.2f\n', IEPSC(kmax(2)), vmax(1), vm(kmax(2)));
%fprintf(1, 'IEPSC(-60) = %6.2f  V(-60)=%6.2f %6.2f\n', IEPSC(kmin(2)), vmin(1), vm(kmin(2)));
warning on MATLAB:divideByZero
warning on MATLAB:nearlySingularMatrix

QueMessage('EPSC_Rect - Storing results');

EPSC_Rect.record = CONTROL(sf).recbeg;
EPSC_Rect.V = vm;
EPSC_Rect.IEPSC = IEPSC;
EPSC_Rect.QEPSC = QEPSC;
EPSC_Rect.Iratio = -IEPSC(kmax(2))/IEPSC(kmin(2));
EPSC_Rect.Gratio = imax/imin;
EPSC_Rect.Vrev = vrev;
EPSC_Rect.WHBlocker = p(2);
EPSC_Rect.WHKd0 = p(3);
EPSC_Rect.WHzdelta = p(4);
EPSC_Rect.WHIndex = p(3)/pA;
EPSC_Rect.WHVr = Vr;
EPSC_Rect.WHRTF = RTF;
EPSC_Rect.T = T;
EPSC_Rect.pA = pA;

CONTROL(sf).EPSC_Rect = EPSC_Rect; % copy the structure over

QueMessage('EPSC_Rect - analysis complete');

%[idiff] = FP_artsupp(idiff, DFILE, sf);

%-----------------------------Prepare for plotting--------------------------------
% for plotting, do baselines/stddev.
%----------------------------- plot if figure is set--------------------------------
if(plot_flag >= 0)
   h = findobj('Tag', 'EPSC_Rect'); % check for pre-existing window
   if(isempty(h)) % if none, make one
      h = figure('Tag', 'EPSC_Rect', 'Name', 'EPSC Rectification', 'NumberTitle', 'off');
      datac('addwindow', 'EPSC_Rect');
   end
   figure(h); % otherwise, select it
   clf; % always clear the window...
   fsize = 7;
   msize = 3;
   tmax = max(time);    
   
   
	set(gca, 'FontSize', fsize);
   ylabel('EPSP Amplitude (nA)', 'FontSize', fsize);
   set(gca, 'XTickLabelMode', 'Manual');
   %grid;
   u=get(gca, 'YLim');
   % plot the EPSPs...
   subplot('Position', [0.1, 0.5, 0.4, 0.4]);
   plot(time([2:2:end],tep1:tep2)', idiff(:,tep1:tep2)'/1000);
   xlabel('ms', 'FontSize', fsize);
   ylabel('nA', 'FontSize', fsize);
	set(gca, 'FontSize', fsize);
   
   subplot('Position', [0.55, 0.5, 0.4, 0.4]);
   plot(vm, IEPSC/1000, 'ko', 'markerfacecolor', 'black', 'markersize', 2);
   hold on;
   plot(vm, Vp/1000, 'r-');
   plot(vm, ywoodhull/1000, 'gs-');
   plot([min(vm) max(vm)], [0 0], 'k-');
   plot([0 0], [min(IEPSC/1000) max(IEPSC/1000)], 'k-');
   xlabel('Vm (mV)', 'Fontsize', fsize);
   ylabel('nA (mV)', 'Fontsize', fsize);
   set(gca, 'Fontsize', fsize);
   

   % plot textual information abstracted from analysis...
   subplot('Position',[0.05,0.05,0.9,0.30])
  	axis([0,1,0,1])
   axis('off')
   h=text(0, 0.9, sprintf('%-12s R[%d:%d] Protocol:%-8s',DFILE.filename, DFILE.frec, DFILE.lrec, CONTROL(sf).protocol), 'Fontsize', 8);
   set(h, 'Interpreter', 'none');
   h=text(0, 0.8, sprintf('Comment: %s', DFILE.comment), 'Fontsize', 8);
   set(h, 'Interpreter', 'none');
   text(0, 0.7, sprintf('Gain:%4.1f LPF:%4.1f kHz Sol:%-12s  ', DFILE.igain, DFILE.low_pass(1), CONTROL(sf).solution), 'Fontsize', 8);
   text(0, 0.6, sprintf('Rectification Ratios: I=%6.4f G=%6.4f (I-60 = %6.2f @ V=%6.2f mV,  I+40 = %6.2f @ V=%6.2f mV)', ...
   EPSC_Rect.Iratio, EPSC_Rect.Gratio, IEPSC(kmin(2)), vmin, IEPSC(kmax(2)), vmax), 'Fontsize', 8); 
   h=text(0, 0.5, sprintf('Reversal Potential: %6.2f mV  (JP = %6.1f mV)  Rs: %7.1f Mohm Comp: %7.1f', ...
       vrev, number_arg(CONTROL(sf).jpot), number_arg(CONTROL(sf).access), number_arg(CONTROL(sf).comp)), 'Fontsize', 8);
   text(0, 0.4, sprintf('Woodhull Parameters: [Blocker] = %8.3f  Kd(0) = %8.3f   z(1-delta) = %8.3f', p(2), p(3), p(4)), 'Fontsize', 8);
    text(0, 0.3, sprintf('Rectification Index (KD(0)/[pa]) = %8.3f  Vr = %8.2f', p(3)/pA, Vr), 'Fontsize', 8);
    text(0, 0.2, sprintf('Temperature: %6.1f   [pa]: %8.1 uM  Fit error: %8.3f', T, pA, e), 'Fontsize', 8);

   set(h, 'Interpreter', 'tex');
   orient landscape
   drawnow
   % control printing and closing of window for automatic runs
   % f = 1 creates plot and leaves it up
   % f = 2 creates plot but closes it when done
   % f = 3 creates plot and prints it and then closes it
   if (plot_flag > 0)
      print -dljet3;
   end
   if plot_flag == 2  
      close
   end
end;

